</div>
Examen 3
Física Computacional
Autor: Celeste Castro Granados $\mathbb C \hat{e}l \mathbb s$
using Plots
using PyPlot
using LaTeXStrings
using LinearAlgebra
using DataFrames
1. Encuentra la solución numérica para la ecuación diferencial $\nabla^2 \phi = 0$ en dos dimensiones con condiciones a la frontera $\phi(0, y) = \phi(x, 0) = \phi(1, y) = \phi(x, 1) = 0$ y $\phi = 1$ para todos los puntos $x$, $y$ en dos círculos de radio $1/4$ con centros en $(0.4, 0.5)$ y $(0.6, 0.5)$.¶
Dibuje las líneas de campo vectorial que se generan de obtener $−\nabla \phi$.¶
a) Solución numérica¶
Para este problema, vamos a utilizar la discretización en 2D, por lo que, comenzaremos programando la función para construir la malla que necesitamos, en la cual $h_x=h_y=h$:
function Malla_cyr_2(x,y,r̃,V₀=1.0) #malla para dos circulos
Mat = zeros(length(y),length(x)) # indices i,j
for i in 2:length(y)-1 #el for va de 2:length(y)-1 porque no estamos tomando la frontera
for j in 2:length(x)-1
#fijamos los centros de ambos círculos
if sqrt( (x[j]-0.4)^2 + (y[i]-0.5)^2 ) <= r̃ || sqrt( (x[j]-0.6)^2 + (y[i]-0.5)^2 ) <= r̃
Mat[i,j] = V₀ #valor para los puntos interiores de ambos círculos
else
Mat[i,j] = rand() #valor para los puntos que no están dentro de los círculos
end
end
end
return Mat
end
Y definimos nuestra malla (que va a tener los valores iniciales):
x = collect(0:0.01:1)
y = collect(0:0.01:1)
mi_malla = Malla_cyr_2(x,y,0.25)
heatmap(x,y,mi_malla,aspect_ratio=:equal)
Ahora, definimos la función que utilizaremos para relajar la malla, la cual utiliza el método de relajación:
function relax!(Mat,V₀=1.0)
for i in 2:length(y)-1
for j in 2:length(x)-1
if Mat[i,j] ≠ V₀
Mat[i,j] = 0.25*( Mat[i+1,j] + Mat[i-1,j] +
Mat[i,j+1] + Mat[i,j-1] )
end
end
end
return Mat
end
Luego, probaremos la función con nuestra malla definida arriba:
mi_malla = relax!(mi_malla)
heatmap(x,y,mi_malla,aspect_ratio=:equal)
Observamos que la función relax! relaja la malla una vez cada que se ejecuta, sin embargo, una vez no es suficiente para que la malla converja a la solución y tendríamos que estar ejecutando dicha función hasta que se alcance la convergencia, por lo tanto, necesitamos otra función que relaje la malla las veces necesarias, y nos entregue como resultado la malla una vez que ya haya convergido. Para esto, ocuparemos el método de Jacobi:
function met_jacobi_v1!(Mat,error=1e-5)
testigo, cuenta = true,0
while testigo == true
Mat_old = copy(Mat)
Mat = relax!(Mat) #va a ir relajando la malla en cada vuelta del while
cuenta += 1
if maximum( abs.( Mat .- Mat_old ) ) <= error #condición para que termine el proceso de relajación
testigo = false
end
end
return Mat, cuenta
end
Finalmente, obtenemos ahora sí la solución:
mi_malla = Malla_cyr_2(x,y,0.25)
mi_malla,cuenta = met_jacobi_v1!(mi_malla)
heatmap(x,y,mi_malla,aspect_ratio=:equal) #la malla ya convergió a la solución
Y para terminar, podemos graficar también las equipotenciales utilizando la función contour:
p1 = contour(x, y, mi_malla,nlevels=10,fill=true)
plot(p1,aspect_ratio=:equal)
b) Líneas de campo vectorial que se generan de obtener $−\nabla \phi$.¶
Para este inciso, necesitamos obtener el campo eléctrico, i.e. el gradiente de la solución del inciso anterior. Por lo tanto, vamos a definir primero una función que calcule justamente el gradiente de una función:
#Definiremos dos funciones, una que regrese la parte x del gradiente y otra que regrese la parte y
#Parte x
function gradiente_x(x,y,ϕ) #ϕ-función a la que se le va a calcular el gradiente
h=x[2]-x[1]
grad_x=zeros(length(x),length(y))
for i in 2:length(y)-1
for j in 2:length(x)-1
grad_x[i,j]=(ϕ[i,j+1]-ϕ[i,j-1])/h #calculo de la derivada
end
end
return grad_x
end
#Parte y
function gradiente_y(x,y,ϕ)
h=y[2]-y[1]
grad_y=zeros(length(x),length(x))
for i in 2:length(y)-1
for j in 2:length(x)-1
grad_y[i,j]=(ϕ[i+1,j]-ϕ[i-1,j])/h #calculo de la derivada
end
end
return grad_y
end
Y con lo anterior, ya podemos graficar el campo eléctrico, para lo cual utilizaremos la función streamplot:
E_x=-1 .*gradiente_x(x,y,mi_malla)
E_y=-1 .*gradiente_y(x,y,mi_malla)
streamplot(x,y,E_x,E_y)
2. Encuentre la solución a la ecuación de Fourier en 2D con condición a la frontera $\phi(0, y) = \phi(x, 0) = \phi(1, y) = \phi(x, 1) = 0$ y con condición inicial $\phi(x, y, 0) = 0$ para todos los puntos afuera de los dos círculos de radio $1/4$ con centros en $(0.4, 0.5)$ y $(0.6, 0.5)$, y $\phi = 1$ para todos los puntos dentro del los dos círculos anteriores.¶
Puede suponer que el coeficiente de difusión es 1.¶
Vamos a definir primero una función que nos permita obtener la matriz de evolución para poder obtener la solución del problema siguiendo el método de Cranck-Nicholson como lo vimos en clase:
function Mat_evol_CN_2D(x,y,t,D=1.0)
Nx,Ny=length(x),length(y)
Δx,Δy,Δt=x[2]-x[1],y[2]-y[1],t[2]-t[1]
r = (D*Δt)/(2*Δx*Δy)
println("r = ",r)
A=zeros( (Nx-2)*(Ny-2) ,(Nx-2)*(Ny-2) )
B=zeros( (Nx-2)*(Ny-2) ,(Nx-2)*(Ny-2) )
#Para la diagonal
for i in 1:Nx-2
for j in 1:Ny-2
A[ (Ny-2)*(i-1)+j , (Ny-2)*(i-1)+j ]=1+4*r
B[ (Ny-2)*(i-1)+j , (Ny-2)*(i-1)+j ]=1-4*r
end
end
#Para las diagonales inferior y superior
for i in 1:Nx-2
for j in 1:Ny-3
A[ (Ny-2)*(i-1)+j , (Ny-2)*(i-1)+j+1 ]=-r
A[ (Ny-2)*(i-1)+j+1 , (Ny-2)*(i-1)+j ]=-r
B[ (Ny-2)*(i-1)+j , (Ny-2)*(i-1)+j+1 ]=r
B[ (Ny-2)*(i-1)+j+1 , (Ny-2)*(i-1)+j ]=r
end
end
for i in 1:Nx-3
for j in 1:Ny-2
A[ (Ny-2)*(i-1)+j , (Ny-2)*(i)+j ]=-r
A[ (Ny-2)*(i)+j , (Ny-2)*(i-1)+j ]=-r
B[ (Ny-2)*(i-1)+j , (Ny-2)*(i)+j ]=r
B[ (Ny-2)*(i)+j , (Ny-2)*(i-1)+j ]=r
end
end
A_inversa=inv(A)
Evol = A_inversa*B
return Evol
end
Y utilizando la función anterior definimos el método de Cranck-Nicholson para 2 dimensiones:
function Cranck_Nicholson_2D(Ψ₀,M_evo,t)
#Vamos a guardar las diagonales
Nx,Ny,Nt=length(Ψ₀[:,1]),length(Ψ₀[1,:]),length(t)
Ψ_xy_t=zeros(Nx,Ny,Nt) #Copiando Psi inicial a la solución final
Ψ_xy_t[:,:,1] = Ψ₀
Ψ_in_t = reshape( Ψ₀[2:end-1,2:end-1],((Nx-2)*(Ny-2),1)) #Vector
for i in 2:Nt
Ψ_in_t = M_evo*Ψ_in_t
Ψ_xy_t[2:end-1,2:end-1,i]=reshape(Ψ_in_t,(Nx-2,Ny-2))
end
return Ψ_xy_t
end
Luego, vamos a definir otra función más, que será la que va a construir a la $\Psi_{inicial}$:
function Psi_ini_circ(x,y)
Nx,Ny = length(x),length(y)
phi=zeros(Ny,Nx)
#cordenada y
for i in 2:Nx-1
for j in 2:Ny-1
if sqrt(( x[i]-0.5 )^2+(y[j]-0.5)^2)<=(1/4)
phi[j,i]=1
end
end
end
return phi
end
Ahora sí, ya tenemos todo lo necesario para resolver el problema. Definimos primero la $\Psi_{inicial}$ a partir de la función Psi_ini_circ:
x=collect(0:0.02:1)
y=collect(0:0.02:1)
t=collect(0:0.0001:0.05)
psi_mat = Psi_ini_circ(x,y)
heatmap(x,y,psi_mat,aspect_ratio=:equal) #valores iniciales
Luego, definimos la matriz de evolución con la función que programamos arriba y obtenemos la solución ($\Psi_{total}$) siguiendo el método de Cranck-Nicholson:
M_t= Mat_evol_CN_2D(x,y,t)
Psi_total = Cranck_Nicholson_2D(psi_mat,M_t,t)
heatmap(x,y,Psi_total[:,:,200],clim=(0,1),aspect_ratio=:equal) #valores finales
Y calculamos las isotermas:
contourf!(x,y,Psi_total[:,:,200],clim=(0,1),aspect_ratio=:equal)
Las dos últimas gráficas representan ya la solución con los valores finales, por lo tanto, vamos a generar también una animación que muestre como va evolucionando la matriz conforme transcurre el tiempo, hasta llegar a la solución final, para poder observar todo el proceso (proceso de difusión del calor):
mi_peli=@animate for i =1:length(t)
heatmap(x,y,Psi_total[:,:,i],clim=(0,1),aspect_ratio=:equal)
end
gif(mi_peli,"2circ_calor.gif",fps=20)
a) Estados estacionarios¶
Vimos en clase que podemos ver este problema como un problema de eigenvalores matricial, por lo que vamos a definir primero una función que nos construya la matriz que necesitamos:
function Mat_onda_2D(x,y)
Nx,Ny = length(x),length(y)
A = zeros(Ny*Nx,Ny*Nx)
for i in 1:Nx
for j in 1:Ny
A[Ny*(i-1)+j,Ny*(i-1)+j] = -4
end
end
for i in 1:Nx
for j in 1:Ny-1
A[Ny*(i-1)+j,Ny*(i-1)+j+1] = 1
A[Ny*(i-1)+j+1,Ny*(i-1)+j] = 1
end
end
for i in 1:Nx-1
for j in 1:Ny
A[Ny*(i-1)+j,Ny*i+j] = 1
A[Ny*i+j,Ny*(i-1)+j] = 1
end
end
return A
end
Ahora, con esta función definimos nuestra matriz, considerando que $L_x = \pi$ y $L_y = 2\pi$:
x = collect(0:0.05:π) #Lx=π, h=0.05
y = collect(0:0.05:2*π) #ly= 2π, h=0.05
M = Mat_onda_2D(x[2:end-1],y[2:end-1]); #la matriz va a estar construida solo con los puntos interiores
Y vamos a obtener los eigenvalores de la matriz M utilizando la función eigvals.
e_vals_M = eigvals(M);
Luego, siguiendo el desarrollo del problema que hicimos en clase, observamos que los valores que corresponden a los estados estacionarios, que son los eigenvalores que estamos buscando realmente, están dados por:
$$e_{vals} = \frac{\sqrt{-\alpha}}{h} $$En donde, $h$ es la separación en la malla (en la cual estamos suponiendo que $h_x=h_y=h$) y las $\alpha$ son los eigenvalores de la matriz M.
h = x[2]-x[1] #hx=hy=h
e_vals=sqrt.(-e_vals_M)/(h);
Finalmente, vamos a comparar los valores que acabamos de obtener (eigenvalores numéricos) con los eigenvalores esperados analíticamente. Para esto, recordamos que los eigenvalores para este problema están dados por:
$$ e_{vals}= \pi \sqrt{\frac{n^2}{L_x^2} + \frac{m^2}{L_y^2}} \leftarrow \text{Valores analíticos}$$En donde, $n$ y $m$ son números enteros.
Lx=π
Ly=2*π
e_vals_a=[]
#Vamos a comparar solo los primeros 25 eigenvalores
for n in 1:5
for m in 1:5
e_val = π * ( sqrt.( (n^2/Lx^2) + (m^2/Ly^2) ) )
push!(e_vals_a,e_val)
end
end
e_vals_n=sort(e_vals) #Odenamos de menor a mayor el arreglo de los eigenvalores
DataFrame(e_valores_analíticos=e_vals_a,e_valores_numéricos=e_vals_n[1:length(e_vals_a)])
Por lo tanto, los eigenvalores encontrados numéricamente si coinciden con los analíticos al menos en la primera cifra (en la mayoría). Para alcanzar una mayor precisión tendríamos que disminuir el paso en los arreglos de $x$ y $y$, sin embargo, con un paso de $0.01$ se detuvo el proceso por falta de memoria, por lo que se decidió dejar el paso en $0.05$.
b) $L_x = \pi$, $L_y =\sqrt{2\pi}$¶
Repetimos el procedimiento anterior, pero ahora para estos valores de $L_x$ y $L_y$:
x = collect(0:0.05:π)
y = collect(0:0.05:sqrt(2*π))
M = Mat_onda_2D(x[2:end-1],y[2:end-1])
e_vals_M = eigvals(M)
h = x[2]-x[1] #hx=hy=h
e_vals=sqrt.(-e_vals_M)/(h);
Y nuevamente comparamos con los eigenvalores esperados analíticamente:
Lx=π
Ly=sqrt(2*π)
e_vals_a=[]
#Vamos a comparar solo los primeros 25 eigenvalores
for n in 1:5
for m in 1:5
e_val = π * ( sqrt.( (n^2/Lx^2) + (m^2/Ly^2) ) )
push!(e_vals_a,e_val)
end
end
e_vals_n=sort(e_vals) #Odenamos de menor a mayor el arreglo de los eigenvalores
DataFrame(e_valores_analíticos=e_vals_a,e_valores_numéricos=e_vals_n[1:length(e_vals_a)])
Para este caso, los eigenvalores numéricos también coinciden con los analíticos, aunque solo en la primer cifra (en la mayoría), i.e. con poca precisión, igual que en el inciso anterior.